
//****************************************************************************************
// GM_ERBSeval.h: interface for the ERBSeval class.
//****************************************************************************************
 
#ifndef GM_ERBSEVAL_H_
#define GM_ERBSEVAL_H_

#include <cmath>
#include <algorithm>
 
class ERBSeval
{
      struct Double5
      {
            double p[5];
            double& operator[](int i) {return p[i];} 
      };
 
      int               j;
      double            dt;
      double            _dt, _sk, _sk1, _sk2;
      int               _m, _der;
      double            _tk, _dtk;
 
      double*           _b;
      Double5*					_a;
 
      double            alpha, beta, gamma, lambda;
 
 
      double phi(double t)
      {
            double d = pow(t*pow(1-t,gamma), alpha);
            if(d < 2.3e-308)
                  return 0.0;
            else 
                  return exp(-beta*pow( fabs(t-lambda), alpha*(1+gamma))/d);
      }
 
      double f2(double t)
      {
            if(t==lambda || t<2.3e-308 || t==1.0)    return 0.0;
            
            double h = (t-1/(1+gamma))*fabs(t-lambda)/(t*(1-t));
            if(t<lambda)      h -= 1.0;
            else              h += 1.0;
            h *= -_sk*alpha*beta*(1+gamma)/pow(t*pow(1-t,gamma),alpha);
 
            if((1+gamma)*alpha < 1)
            {
                  double g = pow(fabs(t-lambda),1-(1+gamma)*alpha)/h;
                  if(g < 2.3e-308)  return 0.0;
                  else                    return 1/g;
            }
            else if ((1+gamma)*alpha > 1)
                  return      h * pow(fabs(t-lambda),(1+gamma)*alpha-1);
            else
                  return      h;
      }
 
      double integral(double a, double b, double sum, double eps)
      {
            double M[16][16];
            int i,j;
            double t   = b-a;
            M[0][0]    = sum*t;
            for(i=1; i<16; i++) 
            {
                  double s = 0.0;
                  int    k = 1 << i;
                  t /= 2.0;
                  for (j=1; j<k; j+=2)  s += phi(a+t*j); 
                  M[0][i] = t*(sum += s);
                  for (j=1; j<=i; j++)
                  {
                        b = 1 << (j << 1);
                        M[j][i-j] = (b*M[j-1][i-j+1]-M[j-1][i-j])/(b-1.0);
                  }
                  if(abs(M[i][0]-M[i-1][0]) < eps) return M[i][0];
            }
            return M[15][0];
      }
 
      
      void interpolate(int i, double p0, double p1, double f0, double f1)
      {
            _a[i][0] = p0;
            _a[i][1] = f0;
            _a[i][2] = 3*(p1-p0) - 2*f0 - f1;
            _a[i][3] = -2*(p1-p0) + f0 + f1;
            _a[i][4] = (p0+p1)/2 + (f0-f1)/12;
      }
 
 
      void _init(int m)
      {
            int i,j;
            double p0, p1;
            double f0, f1;
            _b    = new double[m+1];
            _a    = new Double5[m];
            _m    = m;
            _dt   = 1.0/m;
            _sk   = 1.0;
            _b[0] = p0 = f0 = 0.0;
            for(i=1; i<m; i++)
            {
                  double t = _dt*i;
                  p1    = phi(t);
                  f1    = f2(t)*p1;
                  _b[i] = integral(t-_dt, t, 0.5*(p0+p1), 1e-17);
                  interpolate(i-1, p0, p1, _dt*f0, _dt*f1);
                  p0    = p1;
                  f0    = f1;
            }
            _b[m] = integral(_dt*(m-1), 1, 0.5*p0, 1e-17);
            interpolate(m-1, p0, p1, _dt*f0, 0.0);
            for(i=1; i<=m; i<<=1)
                  for(j=m; j>=i; j-=1)
                        _b[j] += _b[j-i];
            _sk = 1/_b[m];
      }
 
 
      void prepare(double t)
      {
            t  = (t-_tk)/(_dtk);                     // translate/scale parameteren 1
						j   = std::min<int>(int(t*_m),_m-1);               // find local interval
            dt  = (t - j*_dt)/_dt;                   // translate/scale parameteren 2
      }
 
public:
      ERBSeval(int n=1024, double a=1, double b=1, double g=1, double l=0.5)
      {
            _der = 1;
            set(a,b,g,l);
            set_t(0,1);
            _init(n);
      }
 
      virtual
      ~ERBSeval()
      {
            delete _b;
            delete _a;
      }
 
      double operator()(double t)
      {
            prepare(t);
            if(dt > 0.5)
							return 1.0 - (_sk*(_b[j+1] - _dt*(_a[j][4]-dt*(_a[j][0]+dt*(_a[j][1]/2+dt*(_a[j][2]/3+dt*_a[j][3]/4)))))); 
            else
							return 1.0 - (_sk*(_b[j]   + _dt*dt*(_a[j][0]+dt*(_a[j][1]/2+dt*(_a[j][2]/3+dt*_a[j][3]/4)))));
      }
 
      double der1(double t)   { prepare(t); return der1(); }
      double der2(double t)   { prepare(t); return der2(); }
 
      double der1()     { return - _sk1*(_a[j][0] + dt*(_a[j][1]+dt*(_a[j][2]+dt*_a[j][3]))); }
      double der2()     { return - _sk2*(_a[j][1] + dt*(2*_a[j][2]+dt*3*_a[j][3])); }
      double Sk()       { return _sk; }
 
      double der(int i) 
      {
            switch(i)
            {
							case 1: return der1();
							case 2: return der2();
            }
      }
 
      void set(double a=1, double b=1, double g=1, double l=0.5)
            { alpha=a; beta=b; gamma=g; lambda=l; }
 
      void set_t(double tk, double dtk) 
            {_tk=tk; _dtk=dtk; _sk1 = _sk/_dtk; _sk2=_sk/(_dt*_dtk*_dtk);}
};
 
#endif // GM_ERBSEVAL_H_
